1/ / Noise and Long Configuration Memory 
in Bak-Tang-Wiesenfeld Models on Narrow Stripes 
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We report our findings of an 1// power spectrum for the total amount of sand in directed and undi- 
rected Bak-Tang-Wiesenfeld models confined on narrow stripes and driven locally. The underlying 
mechanism for the 1// noise in these systems is an exponentially long configuration memory giving 
rise to a very broad distribution of time scales. Both models are solved analytically with the help 
of an operator algebra to explicitly show the appearance of the long configuration memory. 



PACS number(s): 05.65.+b, 05.45.-a, 05.40.-a 
I. INTRODUCTION 

The ubiquitous 1// noise fascinated physicists for gen- 
erations fug] . There are many examples from wildly dif- 
ferent systems in which the power spectra S(f) ~ 1//" 
with a close to one. This phenomenon usually indicates 
the presense of a broad distribution of time scales in the 
system. The common case where a = 1 is particularly 
intriguing, in that it implies a kind of "equal partition" 
of power among every decade of frequency range, i.e. the 
integral 
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is independent of /. One mechanism to generate such 
a distribution of time scales in systems at thermal equi- 
librium is through thermal activation events over a suffi- 
ciently broad and flat distribution of energy barriers pi . 
The "local" power spectrum generated by any single bar- 
rier has a characteristic frequency which decreases expo- 
nentially with increasing barrier height. But the super- 
position of the power spectra from all the barriers gives 
rise to an 1/ / spectrum. This mechanism is often em- 
ployed to explain, for example, the 1// spectrum of low 
frequency voltage fluctuations in semi-conductors S. In 
search for a more general answer, applicable to nonequi- 
librium and dynamic systems, Bak, Tang, and Wiesen- 
feld (BTW) introduced the notion of Self-Organized Crit- 
icality (SOC) ||. In particular, they proposed a simple 
"sandpile" (BTW) model which shows the emergent scale 
free behavior in both space and time. However, the origi- 
nal BTW model did not exhibit the 1// noise B. In this 
paper we report the observation of 1// noise for directed 
and standard (undirected) BTW models confined on nar- 
row stripes (quasi-one-dimensional geometries). In these 
models, sand flows in the long direction, with periodic or 
closed boundary conditions in the other direction. The 
system is driven locally by randomly adding sand to a 
unique set of sites that have the same coordinate along 
the long axis. The total amount of sand in the sandpile 



as a function of time, measured by the number of added 
grains, exhibits a clean 1/ / power spectrum with an ex- 
ponentially small lower cutoff. Surprisingly, the mecha- 
nism for the 1// noise in this athermal nonequilibrium 
model is rather similar to the above mentioned thermal 
mechanism. In our model the local characteristic fre- 
quency also falls off exponentially as a function of some 
parameter, which in this case is the distance from the 
driving point. 

II. DIRECTED MODELS 

Let us first consider the simpler directed model, defined 
as follows. An integer variable z(x, y) is assigned on every 
site (x, y) of a two-dimensional lattice of size L x x L y 
(1 < x < L x , 1 < y < L y ). Throughout the paper, we 
refer to z(x, y) as the number of grains of sand (or height) 
at the site (x,y). The dynamics consists of the following 
steps: 

(i) Add a grain of sand to a randomly selected site in 
the first column, (l,y): z(l, y) — > z(l, y) + 1. 

(ii) If as a result of the process the height z(x,y) ex- 
ceeds a critical value z c = 2, the site topples and 
three grains of sand are redistributed from this site 
to three of its nearest neighbors up, down, and to 
the right, that is 

z(x, y + 1) — ► z(x, y + 1) + 1, 

z(x,y- 1) -> z(x,y- 1) + 1, 

z(x + l,y) -> z(x + l,y) + 1, 

z(x,y) -> z(x,y) -3. 

(iii) Repeat step (ii) until all sites are stable, i.e. 
z(x,y) < 2 everywhere. This chain reaction of up- 
dates is referred to as an avalanche. 

(iv) When the avalanche is over, measure the total 
amount of sand in the system 2(t) = ^ z(x,y). 
Then go to step (i). 



Notice that the flow of sand is directed to the right along 
the x-axis. Also note the separation of time scales - the 
duration of individual avalanches is taken to be much 
faster than the unit time interval which is defined by the 
addtion of sand grains. The boundary condition in the x 
direction is always set to open: z(L x + 1, y) — 0. While in 
the y direction we either choose the periodic boundary 
condition (which we refer to as Model 1) or the closed 
boundary condition. In the latter case we restrict our- 
selves to L y — 2 and refer to it as Model 1A. In Model 1A, 
we set z c = 1 and the redistribution rule (ii) prescribes 
to move two grains of sand from the toppling site: one 
to the right along the x direction and the other to its 
nearest neighbor (up or down) in the y direction. 
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FIG. 1. The total amount of sand in the system Z(t) (dots), 
and its running averages over 10 (dashed line), 100 (solid line), 
and 1000 (thick solid line) time steps. The data was taken 
from Model 1A with L x = 8. 

After some transient period the above dynamics brings 
the system to a stationary state, where the total amount 
of sand in the system Z saturates and fluctuates about 
its average value. At this point we start recording Z(t) 
and measure its power spectrum S(f) — \Z(f)\ 2 , where 
Z(f) is the Fourier transform of Z(t). A time window of 
Z(t) is shown in Fig. |l] together with its running averages. 
Notice the fluctuations on many time scales. In Fig. g, 
we show the power spectra for Models 1 and 1A. Even 
for small systems one observes a very broad 1/ / region. 
In fact, as we will demonstrate later, the lower cutoff of 
the 1// region falls off with L x exponentially. Our simu- 
lations indicate that as the width of the stripe in Model 
1, Ly is increased the 1// region shrinks, as shown in 
Fig. 0(c) . Direct observations of configurational changes 
at each time step clearly indicate that the rate of config- 
urational changes at x decreases drastically with increas- 
ing x H , suggesting that there are many time scales and 
some kind of long memory in the system. To understand 
this, we proceed with solving Model 1A using the group 



of operators introduced by Dhar |3||7]] to treat sandpile 
models. 
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FIG. 2. Power spectra for (a) Model 1A with L x 
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and let L x — L. Any pair configuration with 
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Model 1 with L x = 8 and L y = 4; and (c) Model 1 with 
L x = 8 and L y — 8. The curves in (b) and (c) are shifted 
vertically by -1 and -2 decades, respectively, for clarity. The 
dashed line has the slope — 1. 

To simplify the notation let us denote the configura- 
tion at the pair of sites z(x, 1) and z(x, 2) by the column 

z(x,2)' 

z(x,l) 

both z{x,y) < z c = 1 is stable. However, the recurrent 
pair configurations, present in the stationary SOC state, 

!)' (o)' and (i)' while (o 

ized in the system after some transient period. As usual 
in directed models H, there are no additional restric- 
tions on pair configurations at different columns. The 
total number of recurrent (SOC) states is thus 3 L . Let us 
define operators U x and D x acting on recurrent configu- 
rations. The action of these operators consists of adding 
one grain of sand at sites (x, 2) and (x, 1) correspond- 
ingly, and, if necessary, relaxing the resulting configura- 
tion according to the avalanche rules of the dynamics. 
The final stable configuration is the result of the opera- 
tor acting on the initial configuration. One can demon- 
strate that these operators commute with each other, i.e. 
[U X ,U X >] = [D X ,D X ,] = [U X ,D X >] = 0, and the model is 
therefore an Abelian model ||. The following operator 
identities || result directly from the relaxation rules of 
the model 



V. 



L> x U x +i 



D x = U X D X+1 . 



(1) 
(2) 



These identities simply state that the addition of two 
grains of sand to any site (of a recurrent state) will cer- 
tainly make it unstable and, therefore, two grains will 



be transferred to the neighbors according to the relax- 
ation rule of Model 1A. Open boundary at x = L + 1 
corresponds to Ul+i = -Dl+i = I, where I is the identity 
operator. From Eqs. (P and (g) it immediately follows 
that 



U-lDx = U 2 D 2 



U L D L = U L+1 D L+1 = I. (3) 



In other words, addition of one grain of sand to both 
upper and lower sites at any column x triggers a down- 
stream avalanche, in which two grains fall off the right 
edge, but the underlying configuration of sand columns 
remains unchanged. Since D x = U x l : Eq. (fil) can be 
rewritten as Rjj = U x+ \. Therefore, Uf — U$ = 
U\ = U L+1 , and 



Uf =1. 



(4) 



Repeated application of U\ 3 L times makes the system 
visit every one of the 3 recurrent states exactly once 
and return to its original configuration. If U\ and D\ 
are applied in random order (that is how we drive the 
system in Model 1A), eventually all recurrent states will 
still be visited. But, since [/{"£>" = U™~ n , for a random 
sequence of C/i's and £>i's the average time required to 
visit all 3 L states is given by VT = 3 L , or T = 9 L . Now 
let us study the "microscopic" details of how operators 
U x and D x change the configurations of the sandpilc. 
Note that 
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(5) 
(6) 

(7) 



It is clear that the only cases where the operator U x or 
D x can "propagate" to the right and change the state 

of the next neighbor are U x ( n ) = ( 1 j < . .... ., 



and 
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D., 
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So, in order for the operator 



Thus, the probability that a pair at x changes its 



U\ to be able to change the configuration of a pair at x, 
all x — 1 pairs to the left of this pair must be in the state 



configuration as a result of the addition of a grain of sand 
at the left end of the system is 1/3 X ~ 1 . In other words, 
the system has an exponentially long configuration mem- 
ory. If the system is driven by a random sequence of U\ 
and £>x, on average it will take t x ~ (3 X ~ 1 ) 2 = 9 X ~ 1 
grains of sand to change the configuration of a pair at x. 

This exponential increase of the characteristic time t x 
with x is manifested in the local autocorrelation func- 
tions C(x,t) = [< Z(x,0)Z(x,t) > - < Z(x,0) > 2 ]/[< 
Z{x,Q) 2 > - < Z(x,0) > 2 ], where Z(x,t) = z(x,l) + 



z(x,2) is the number of grains at column x and at time 
t. In Fig. H we plot C(x, t) vs. t/9 x_1 for several x's. One 
sees that C(x, t) = F(t/9 a:_1 ). This form of the local au- 
tocorrelation function implies a scaling form for the local 
power spectrum: Si oc (f,x) = (l/f c har(x))S(f/f char (x)) 
with f c har(x) = / Gxp(-a;ln9), where Si oc {f,x) is the 
power spectrum of Z(x,t). Note that (see Eq. (H)) if U x 
or D x propagates through a column it leaves the num- 
ber of grains on that column Z(x) unchanged. It follows 
that the addition of a grain at the left end of the sys- 
tem can change at most the number of grains at one 
column. If we assume that the local events of changing 
Z[x) are independent for different x's (or the correla- 
tions of which is not too strong), which is a reasonable 
approximation when we drive Model 1A with a random 
sequence of E/i's and -Di's, then the global power spec- 
trum of the total number of grains in the system is the su- 
perposition of the local power spectrum. The exponential 
fall off of the local characteristic frequencies of configura- 
tion changes would give rise to a global 1// power spec- 
trum, similar to the case of thermal activations in equi- 
librium systems mentioned earlier S. That is S g (f) — 

Jo Sl oc (f,x)dx = J (1/ fchar(x))S(f / fchar(x))dx = 

J L eM^)S(feM^)/fo)dx/f = (l/f)fffi lXL)/fo 
dyS(y)/X. The lower cutoff of the 1// region is f c ~ 
/oexp(— XL), which for the curve (a) in Fig. g (A = In 9 
and L = 8) is of the order 10~ 7 . 
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FIG. 3. Autocorrelation functions C(x, t) at X — 3, x — 4, 
and x — 5, for Model 1A with L x = 8. 



III. UNDIRECTED MODELS 

We now turn our attention to the undirected sandpile 
model on a stripe L x x L y . In this model an unstable 
site with z(x, y) > z c — 3 redistributes one grain of sand 
to each of its four neighbors. In our simulations at each 
time step we randomly select a site on the central column 



(x = (L x + 1)/2, or, if L x is an even number we randomly 
select one of the 2L y sites on the two central columns) 
and add one grain of sand to that site. We choose to have 
open boundaries at x — and x = L x + 1. In the first 
version of the model the boundary condition along the y 
direction is periodic. We refer to this version as Model 
2. In addition to Model 2 we consider a simpler model 
which we refer to as Model 2A. Model 2A is defined on 
an L x 2 stripe with closed boundary condition in the y- 
direction. In Model 2A z c — 2 and a site with z(x, y) > 2 
moves one grain of sand to each of its three neighbors. 
The advantage of studying Model 2A is that it has fewer 
recurrent configurations which are easier to classify. In 
Fig. [|, we show the power spectra of the total amount 
of sand in Models 2 and 2A. Similar to the case of the 
directed Model 1, in Model 2 the 1// region shrinks for 
increasing L y . The dynamics of the undirected models is 
apparently more complex than that of the directed ones. 
However, much of the apparent complexity is due to the 
motion of "troughs" g - columns in which all 2 < z c — 1, 
so that avalanches can not propagate beyond them fl9|. 
Let us first understand the trough dynamics. 
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FIG. 4. Power spectra for the center-driven BTW models: 
(a) Model 2A with L x = 16, (b) Model 2 with L x = 17 and 
L y = 8, and (c) Model 2 with L x = 17 and L y = 24. The 
dashed line has the slope — 1. 



A. Trough dynamics 

Let us focus on Model 2 A. Like in Model 1A, one can 
define operators U x and D x . The operator algebra now 
satisfies the following operator relations: 



U x — D x U x +\U x -\ 
D x = U x D x+ iD x _i 



(8) 
(9) 



operator. Let us define the operator O x = U X D X . In 
Model 1A we have shown that O x = I for every x. This 
is not so in the undirected model. However, in this model 
the operators O x form a simple small subgroup of all 
operators in the system. From Eqs. (||) and (0) it fol- 
lows that O x = O x+ iO x -\. Using this operator identity 
repeatedly one gets O x = O x+ iO x -i = O x+2 O x - 2 = 
. . . = O x + n O x - n . Similarly, one can show that O x O x > = 
O x -iO x >+i = O x -20 x >+2 = . . . = OnO x + x i = O x + x >, for 
x + x' < L. In general, 

O x O x ' = 0( ffi+x /) mod (i+i). (10) 

In other words, operators Oo(= I), 0\, 2 , ■ ■ ■ , Ol form 
a cyclic subgroup of L + 1 elements |l(| . 

To understand the physical nature of this subgroup let 
us take a closer look at the set of recurrent configurations 
in Model 2A. In a stable configuration each z(x,y) can 
take values 0, 1, 2. Each pair thus has 9 stable configura- 
tions. However, the number of recurrent configurations 
is much smaller than 9 L . To check which stable configu- 
rations are forbidden in the recurrent set one applies the 
rule developed in o. According to this rule, a subcon- 
figuration at a subset of sites F is forbidden if for every 
site (x, y) £ F z(x, y) is strictly smaller than the number 
of its neighbors in the subset F . It is clear that the pair 

I is forbidden. Let us refer to pairs I 1 ] , I 

1 



propagating beyond them. It is easy to see that a sub- 
configuration enclosed by two troughs is forbidden. Thus 
a recurrent configuration cannot contain more than one 
trough. Therefore, all SOC states naturally fall into one 
of the L + 1 classes: those with no troughs, and those 
with a trough in the 771-th column. The operator O m 
acting on a state \S m > with a trough at the position m 
does not produce any topplings but simply fills up the 
trough. On the other hand, acting on this state with Ok 
(k =/= m) creates a (usually large) avalanche in which two 
grains of sand fall of the pile. However, this avalanche 
produces only minor changes in the configuration of the 



is forbidden. Let us refer to pairs 

as troughs. Troughs prevent any avalanche from 



pile. Indeed, since Ok = O m O 



-1 



(m-k) mod (L+l) 



the opera- 



tor Ok acting on the state \S m > fills up the trough at 
x = m and creates a trough at X = (m — k) mod (L + 1) 
PJ. One may simply view that Ok moves the trough 
from x = mtox = m~ k mod (L + 1). The action of Ok 
on a state \S > with no troughs results in a system- wide 
avalanche with 4 grains of sand falling off the sandpile. 
The only configurational change, however, is the creation 
of a new trough at L + 1 — k (recall the operator identity 
Ok = Oj] 1 _ k ). These rules mean that the action of the 
L operators Ok results only in the motion (or creation 
and annihilation) of the trough, and does not destroy 
the configuration memory of the system. 



The open boundaries at two ends imply Uq — -Do — 
Ul+i = Dl+i = I, where, as before, I is the identity 



B. The number of recurrent configurations 

Having understood that operators Ok = UkDk, do 
not destroy the long term configuration memory of the 
system but only move, create, and annihilate troughs, 
we proceed with describing how individual operators 
Uu = D^ 1 Ok change the configuration. We separate 
the trivial trough dynamics from others by defining the 
equivalence relation of operators: if A = BOk, we say 
that A is equivalent to B and denote it by A = B. Thus 
Uk — DZ • One can rewrite the basic operator identity 
(§asUl = U k ^U k+1 O k ,or 



U k = U k -iU k+ \. 



(11) 



We now derive the relation between operators U k at dif- 
ferent sites. Since Uq = I, Eq. ( pd| ) implies Uf = U-i- 
Write U[" 1 



N(k) 



u,u. 



kVk+2 



u 



= U k+ \. One has U 1 

N(k-l)jjN(k+l) 



iN(k) 



— u k+l 



, which gives the recursion 
relation N(k + 1) = 4N(k) - N(k - 1) with initial con- 
ditions N(0) = 1, 7V(1) = 4. It is easy to show that 

N(k) = [(3 + 2V3)(2 + V3) fe - (2V3 - 3) (2 - V3) fe ]/6. 

(12) 

In a system of size L one has U x ( ' = Ul+i = I- This 
observation enables us to calculate the total number of 
recurrent (SOC) states as N^ = (L + l)N(L), or 



N, 



(2A) 
SOC 



L+l 
6 



[(3 + 2%/3)(2 + V3) L 
- (2V3-3)(2-V3) L ]. 



(13) 



In other words, any recurrent configuration can be ob- 
tained from a given one by the action of some power 
of U\ (there are N(L) choices of this power before con- 
figurations start repeating themselves) and, if necessary, 
creation, annihilation, or change of the position of the 
trough achieved by the the action of L operators O k . 
Q. Asymptotically, only 2 + \/3 ~ 3.732 pair configura- 
tions per site are allowed in a recurrent state, compared 
to 9 stable pair configurations. The above operator re- 
lations can be easily modified for Model 2 on an L X 2 
lattice. For this model the subgroup ( |l0| ) of operators O k 
remains unchanged, while (11) becomes U® = U k -\U k+ \. 
The number of recurrent configurations in this model is 
given by N ( J% C = (L + 1)[(4 + 3\/2)(3 + 2^2) L - (3^2- 
4)(3-2^) L ]/8. 



C. Configuration memory 

Now we are in the position to address the question of 
long memory in Model 2A. (For Model 2 on an Lx 2 stripe 
these arguments can be repeated step by step with some 
minor changes.) Let us restrict ourselves to the operator 



t/fc acting on a state that has no trough. Indeed, since 
we are interested in general properties of the equivalency 
class produced by the action of O k 's, one can alway select 
from this class a representative state that has no trough. 
The action of the operator U k on the pair at column k is 
given by 



Uk 

u k 

Uk 
Uk 

Uk 



C/fe+it/fe-i 



Oi 



Oi 



(14) 
(15) 
(16) 
(17) 
(18) 



It seems that Eq. ( |l6J ) could propagate U k through a 

string of n I 's, changing the configurations away from 

the driving point, similar to the case of directed mod- 
els. This is not so, because one can not have consec- 
' 2 ' 

x o 

it is Eqs. ([17]) and ([18]) which can cause configuration 
changes away from the driving point. Naively, accord- 



utive columns of 



in a recurrent state. In fact 



ing to Eqs. ([17]) and (|l8|), the action of Uk on 



causes only local changes apart from some 



trough dynamics. However, this is true only if the lo- 
cal changes (l J — > I I or I J — > I I ) do not re- 
sult in any FSC. If an FSC does appear as a result of 
this change, the change in the original configuration will 
not be restricted to one pair, but instead will propagate 
throughout the FSC. It is an easy task to classify all 
FSC's in Model 2A. We have already shown that the pair 

n I and two troughs together with the region between 

them are FSC's. The third important FSC is a string 
(Oil ■••110) (or any other string of 1 's connect- 
ing two 0's, which we disregard for a moment, since it 

must contain at least one trough I , and we currently 

restrict ourselves to the states without troughs.) It is 
easy to check that in such a string z(x,y) everywhere 
is smaller than the number of its neighbors in the FSC. 
Such configuration can be created by the action of Uk on 



Usually one has Uk 



d; 



However, if the 

2 



happens to be in the string of 
s bounded between two ( „ )'s, simply changing 



to 



would create a FSC. Direct application 



2) '" \. I 

of the rules of avalanche dynamics shows that in this 
case all pairs associated with this FSC and the pairs 
next to it will be updated. The other way of creating 
the above mentioned FSC is for Uk to act on a string of 

pairs 0I f x J ended with f Q J , K ^ ^ ^ Q 

Again, in this case all pairs associated with this FSC 
and the pairs next to it will be changed. Both of the 
scenarios require that the starting configuration con- 
tains a string of _. j 's with _ 
the ends. Such a string of length x is just one among 



at least one of 



N. 



(2A) 

soc 



(x) ~ (2 + \f?>) x ~ 3.732 a: recurrent states of a 
string of x columns. That is why the irrevesible changes 
of pairs at distance x from the driving pair are exponen- 
tially unlikely. In Model 2A, driven by random addition 
of sand at sites on the central pair(s) the characteristic 
frequency at distance k from the driving point is given 
by fchar(k) ~ l/((2 + V3) k ) 2 = 13.93~ fc . This analytical 
result is in good agreement with numerical simulations 
of the model shown in Fig. |5|. 
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FIG. 5. Autocorrelation functions C(x,t) for Model 2A 
with L x = 12, at x = 8,9,10,11, which are k = 1,2,3,4 
distance away from the driving pairs x — 6 and x — 7. 



IV. CONCLUSION 

In spite of the apparent differences between directed 
and undirected models, the mechanism for a long term 
memory and the 1/ / spectrum in these two systems is 
the same. We summarize our analysis as follows: 

1) Operators O x = U X D X do not produce irreversible 
changes in the configuration. In the directed mod- 
els these operators do not change the configuration 
at all, while in the undirected models their ability 



to change the configuration is restricted to creation, 
annihilation and motion of the trough. 

2) In order to produce irreversible changes in a con- 
figuration at a distance k from the place of sand 
addition, all k pairs in between have to be in a 
unique peculiar configuration. Since this particular 
subconfiguration is just one among Nsoc(k) ~ A k 
possible recurrent subconfigurations, the character- 
istic frequency of irreversible updates falls off with 
distance exponentially 

3) Such exponential dependence of the characteristic 
frequency of updates leads to the 1/f spectrum of 
the total amount of sand in the sandpile. 

It is straightforward to generate the case to higher dimen- 
sions in which sand flows in one (say the x) direction with 
closed or periodic boundaries in other directions. One 
would still observe the 1/f spectrum. Recently, De Los 
Rios and Zhang |l3|l observed an 1/f spectrum in a non- 
conserved sandpile-like model in which certain fraction of 
sand is lost in each toppling process. Due to the absence 
of conservation avalanches themselves are exponentially 
unlikely to reach a distant site, giving rise naturally to 
an exponential distribution of time scales. In contrast, in 
our model avalanches constantly pass through the system 
but they produce only small changes of the configuration. 
An 1/f spectrum was also observed previously for a con- 
tinuously boundary-driven BTW model [|14| . Its origin 
was attributed to a (linear) diffusion of z(x, y) with a 
noisy boundary condition p3 , 15] , which gives a power- 
law lower cutoff f c ~ 1/Ljfor the 1/f spectrum - a 
mechanism very different than ours. 
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